#!/usr/bin/env R
# -*- coding: utf-8  -*-
#
# This script is free software: you can
# redistribute it and/or modify it under the terms of the GNU General Public
# License as published by the Free Software Foundation, version 2.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc., 51
# Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Copyright Izaskun Mallona
# izaskun.mallona@gmail.com

#
#This parses a XML file produced by MX3000P RealTimePCR
#(chart data, XML horizontally grouped by plot)
#It calculates several efficiency estimates as well CT (Cq, take-off) values
#
#CQ_efficiency and cpCQ_take off are exactly the same as produced
#by the Rotorgene software

library(XML)
library(qpcR)

print("", quote=F)
print("This script requires the XML (horizontally grouped by plot, without ROX) file produced by the MX3000P RealTimePCR machine. It produces a spreadsheet with the crossing points (aka CT, CP or cycle thresholds) and the efficiencies as calculated by several methods:", quote=F)
print("", quote=F)
print("* cpD2, the second derivative method.", quote=F)
print("* cpD1, the first derivative method.", quote=F)
print("* CQ, Corbett Rotorgene's software (it is a cpD2 modified).", quote=F)
print("* Cy0, as described by Guescini M et al, BMC Bioinformatics, 2008, 9: 326.", quote=F)
print("", quote=F)
                        
print("Use example:", quote=F)
print("", quote=F)
print("     eficiencia(fileHandle='C:/R/Chart Data.xml', cycle_number=45)", quote=F)
print("", quote=F)
print("Another use example (it shortens the command):", quote=F)
print("", quote=F)
print("     eficiencia('Chart Data.xml', 45)", quote=F)
print("", quote=F)
print ('Note: Mind the quotation marks or apostrophes. They should flank the input file name.', quote=F)
print ('Note: Mind the path to your file. "C:/R/Chart Data.xml" differs from "C:/tmp/data/Chart Data.xml". Note that it is case sensitive.', quote=F)
print("", quote=F)
print("Finally, REST should receive as inputs CPs and efficiencies obtained with the same method. Cy0 efficiencies pairs cpD2 crossing points.", quote=F)
print("", quote=F)
print("Questions? Complaints? izaskun.mallona@gmail.com", quote=F)
print("", quote=F)

handlerFactory <- function() {
    ## 'local' store visible to functions defined inside
    ## handlerFactory
    counts <- new.env(parent=emptyenv())
    ## return value -- list of functions
    list(startElement=function(name, atts, ...) {
        ## lexical scope often requires use of <<- rather than <-
        if (!exists(name, counts))
            counts[[name]] <- 1
        else
            counts[[name]] <- counts[[name]] + 1
    }, getCounts=function() {
        ## for retrieving results
        as.list(counts)
    })
}

eficiencia<-function(fileHandle, cycle_number) {

	
	sample_number=xmlEventParse(fileHandle,handlerFactory())$getCounts()$sample_master
  
    xml<-xmlTreeParse(fileHandle,useInternal=T)

    children=xmlChildren(xml)[[2]]

    headerReplicate = vector()
    headerName=vector()
    headerType=vector()
    j=8
    dataframe=data.frame()
    data=list()
    y_values=vector()
    while (j < sample_number+8) {
        i=8
        headerReplicate= c(headerReplicate, xmlValue(children[[1]][[j]][[3]]))
        headerName= c(headerName, xmlValue(children[[1]][[j]][[4]]))
        headerType= c(headerType, xmlValue(children[[1]][[j]][[5]]))
        
        while (i < cycle_number+8) {
            y_values = c(y_values, as.numeric(xmlValue(children[[1]][[j]][[7]][[5]][[i]][[6]])))
            
            dataframe= rbind(headerReplicate,headerName,headerType)
            data$id[[j-7]]=dataframe[,j-7]
      
            i=i+1
            
        }
        data$values[[j-7]]=y_values
        merged=c(data$id,data$values)
        merged=cbind(merged,merged)
        y_values=vector()
        j=j+1
    }

    result=list()
    k=1
    while (k <= sample_number) {
        result[[k]]=c(data$id[[k]],data$values[[k]])
        k=k+1
    }

    data=as.data.frame.list(result)

    colnames(data)<-c(1:length(data))

    headerForEfficiency=data[1:3,]
    rownames(headerForEfficiency)= c("Replicate", "Well name", "Well type")    
    data2eff<-apply(as.data.frame(data[4:nrow(data),]),2,as.numeric)
    Cycles<-1:cycle_number
    data2eff<-cbind(Cycles,data2eff)

    i=1
    foo = list()
    while (i < ncol(data2eff)){
        foo[[i]]= c(
            try(efficiency(pcrfit(data2eff,1, i+1),plot=F,type="cpD2")$eff,silent=T),
            try(efficiency(pcrfit(data2eff,1, i+1),plot=F,type="cpD1")$eff,silent=T),
            try(efficiency(pcrfit(data2eff,1, i+1),plot=F,type="CQ")$eff,silent=T),
            try(efficiency(pcrfit(data2eff,1, i+1),plot=F,type="Cy0")$eff,silent=T),
            try(efficiency(pcrfit(data2eff,1, i+1),plot=F,type="cpD2")$cpD2,silent=T),
            try(efficiency(pcrfit(data2eff,1, i+1),plot=F,type="CQ")$cpCQ,silent=T)
        )
                 
        i=i+1 
    }    
    foo2= list()
    k=1
    while (k <=length(foo)) {
		if (length(foo[[k]]) == 6) {
        foo2[[k]]= c(data$id[[k]],data$values[[k]],foo[[k]])
		}
		else
		{
		foo2[[k]] =c(data$id[[k]], data$values[[k]], rep("NA",6))
		}
		
        k=k+1
    }
    
    dataEfficiency=as.data.frame.list(foo2)
    colnames(dataEfficiency)<-c(1:length(data))
    
    result=cbind(t(headerForEfficiency),t(dataEfficiency))
    colnames(result)= c("Replicate", "Well name", "Well type", "cpD2_efficiency", "cpD1_efficiency", "CQ_efficiency", "Cy0_efficiency", "cpD2_takeoff", "cpCQ_takeoff" )
    resultHandle= paste(fileHandle,"statistics.xls")   
    write.table(result,resultHandle,sep="\t", row.names=F)
	
	print("Success. A file called")
	print(resultHandle)
	print("has been created. It's plain text, but you can open it using Microsoft Office, OpenOffice.org Calc or equivalent.")
    }


